\chapter{TOUGH2 thermodynamics}
\label{t2thermo}
\index{thermodynamics!IFC-67}
\index{PyTOUGH!thermodynamics!IFC-67}

\section{Introduction}
The \texttt{t2thermo} library in PyTOUGH contains a Python implementation of the thermodynamic routines used in TOUGH2.  These can be used to calculate the thermodynamic properties of water and steam under a range of conditions.  They are based on a subset of the IFC-67 thermodynamic formulation \citep{IFC_67}.

The \texttt{t2thermo} library can be imported using the command:

\begin{lstlisting} 
   from t2thermo import *
\end{lstlisting}

The functions available through the \texttt{t2thermo} library are listed in Table \ref{tb:t2thermo_functions} and described below.

\index{thermodynamics!IFC-67!functions}
\begin{table}
  \begin{center}
    \begin{tabular}{|l|l|p{65mm}|}
      \hline
      \textbf{Function} & \textbf{Type} & \textbf{Description}\\
      \hline
      \hyperref[sec:t2thermo:cowat]{\texttt{cowat}} & tuple & density and internal energy of liquid water\\
      \hyperref[sec:t2thermo:sat]{\texttt{sat}} & float & saturation pressure as a function of temperature\\
      \hyperref[sec:t2thermo:region]{\texttt{region}} & integer & thermodynamic region\\
      \hyperref[sec:t2thermo:separated_steam_fraction]{\texttt{separated\_steam\_fraction}} & float & separated steam fraction for given enthalpy and separator pressure\\
      \hyperref[sec:t2thermo:supst]{\texttt{supst}} & tuple & density and internal energy of dry steam\\
      \hyperref[sec:t2thermo:tsat]{\texttt{tsat}} & float & saturation temperature as a function of pressure\\
      \hyperref[sec:t2thermo:visw]{\texttt{visw}} & float & dynamic viscosity of water\\
      \hyperref[sec:t2thermo:viss]{\texttt{viss}} & float & dynamic viscosity of steam\\
      \hline
    \end{tabular}
    \caption{\texttt{t2thermo} functions}
    \label{tb:t2thermo_functions}
  \end{center}
\end{table}

\section{Thermodynamic functions}

The thermodynamic routines used in TOUGH2 provide functions for liquid water and dry steam.  These functions calculate secondary parameters from the primary thermodynamic variables.  Their names follow the subroutine names used in the TOUGH2 code.

\begin{snugshade}
\subsection{Liquid water: \texttt{cowat(\emph{t}, \emph{p}, \emph{bounds} = False)}}
\end{snugshade}
\label{sec:t2thermo:cowat}
\index{thermodynamics!IFC-67!liquid water properties}

The \texttt{cowat} function returns a two-element tuple (\texttt{d},\texttt{u}) of density (kg/m$^3$) and internal energy (J/kg) of liquid water as a function of temperature \texttt{t} ($\degree$C) and pressure \texttt{p} (Pa).

\textbf{Parameters:}
\begin{itemize}
\item \textbf{t}: float\\
  Temperature ($\degree$C)
\item \textbf{p}: float\\
  Pressure (Pa)
\item \textbf{bounds}: Boolean\\
  If \texttt{True}, return \texttt{None} if the input temperature and pressure are outside the operating range of the routine (as defined by thermodynamic region 1 of the IFC-67 specification).
\end{itemize}

\begin{snugshade}
\subsection{Dry steam: \texttt{supst(\emph{t}, \emph{p}, \emph{bounds} = False)}}
\end{snugshade}
\label{sec:t2thermo:supst}
\index{thermodynamics!IFC-67!dry steam properties}

The \texttt{supst} function returns a two-element tuple (\texttt{d},\texttt{u}) of density (kg/m$^3$) and internal energy (J/kg) of dry steam as a function of temperature \texttt{t} ($\degree$C) and pressure \texttt{p} (Pa).

\textbf{Parameters:}
\begin{itemize}
\item \textbf{t}: float\\
  Temperature ($\degree$C)
\item \textbf{p}: float\\
  Pressure (Pa)
\item \textbf{bounds}: Boolean\\
  If \texttt{True}, return \texttt{None} if the input temperature and pressure are outside the operating range of the routine (as defined by thermodynamic region 2 of the IFC-67 specification).
\end{itemize}

\section{Viscosity}

\begin{snugshade}
\subsection{Liquid water: \texttt{visw(\emph{t,p,ps})}}
\end{snugshade}
\label{sec:t2thermo:visw}
\index{thermodynamics!IFC-67!viscosity}

The \texttt{visw} function returns the dynamic viscosity (Pa.s) of liquid water as a function of temperature \texttt{t} ($\degree$C), pressure (Pa) and saturation pressure (Pa).

\textbf{Parameters:}
\begin{itemize}
\item \textbf{t}: float\\
  Temperature ($\degree$C)
\item \textbf{p}: float\\
  Pressure (Pa)
\item \textbf{ps}: float\\
  Saturation pressure (Pa), calculated for example using the \texttt{sat} function.
\end{itemize}

\begin{snugshade}
\subsection{Dry steam: \texttt{viss(\emph{t,d})}}
\end{snugshade}
\label{sec:t2thermo:viss}
\index{thermodynamics!IFC-67!viscosity}

The \texttt{viss} function returns the dynamic viscosity (Pa.s) of dry steam as a function of temperature \texttt{t} ($\degree$C) and density \texttt{d} (kg/m$^3$).

\textbf{Parameters:}
\begin{itemize}
\item \textbf{t}: float\\
  Temperature ($\degree$C)
\item \textbf{d}: float\\
  Density (kg/m$^3$)
\end{itemize}

\section{Saturation line: \texttt{sat(\emph{t})} and \texttt{tsat(\emph{p})}}

\begin{snugshade}
\subsection{\texttt{sat(\emph{t}, \emph{bounds} = False)}}
\end{snugshade}
\label{sec:t2thermo:sat}
\index{thermodynamics!IFC-67!saturation curve}

The \texttt{sat} function returns the saturation pressure (Pa) at a given temperature \texttt{t} ($\degree$C), for temperatures below the critical temperature.

\textbf{Parameters:}
\begin{itemize}
\item \textbf{t}: float\\
  Temperature ($\degree$C)
\item \textbf{bounds}: Boolean\\
  If \texttt{True}, return \texttt{None} if the input temperature is outside the operating range of the routine (i.e. less than 0.01$\degree$C or greater than the critical temperature, 374.15$\degree$C ).
\end{itemize}

\begin{snugshade}
\subsection{\texttt{tsat(\emph{p}, \emph{bounds} = False)}}
\end{snugshade}
\label{sec:t2thermo:tsat}
\index{thermodynamics!IFC-67!saturation curve}

The \texttt{tsat} function returns the saturation temperature ($\degree$C) at a given pressure \texttt{p} (Pa), for pressures below the critical pressure.

Note that the IFC-67 formulation did not include an explicit formula for calculating saturation temperature as a function of pressure, so here (as in TOUGH2) this is calculated using an iterative root-finding process on the \texttt{sat} function.  The root-finding function is from the \texttt{scipy} library, so this library must be installed before the \texttt{tsat} function will work.

\textbf{Parameters:}
\begin{itemize}
\item \textbf{p}: float\\
  Pressure (Pa)
\item \textbf{bounds}: Boolean\\
  If \texttt{True}, return \texttt{None} if the input pressure is outside the operating range of the routine (i.e. less than \texttt{sat(0.01)} or greater than the critical pressure, 22.12 MPa).
\end{itemize}

\section{Other functions}

\subsection{Separated steam fraction}

\begin{snugshade}
\subsubsection{\texttt{separated\_steam\_fraction(\emph{h}, \emph{separator\_pressure}, \emph{separator\_pressure2} = None)}}
\end{snugshade}
\label{sec:t2thermo:separated_steam_fraction}
\index{thermodynamics!IFC-67!separated steam fraction}

Returns the separated steam fraction for a given enthalpy \texttt{h} and separator pressure.  A second separator pressure may be specified in the case of two-stage flash.

\textbf{Parameters:}
\begin{itemize}
\item \textbf{h}: float\\
  Enthalpy (J/kg)
\item \textbf{separator\_pressure}: float\\
  Steam separator pressure (Pa)
\item \textbf{separator\_pressure2}: float (or \texttt{None})\\
  Second separator pressure (Pa) for two-stage flash -- set to \texttt{None} (the default) for single-stage.
\end{itemize}

\subsection{Determining thermodynamic region}

\begin{snugshade}
\subsubsection{\texttt{region(\emph{t}, \emph{p})}}
\end{snugshade}
\label{sec:t2thermo:region}
\index{thermodynamics!IFC-67!region}

Returns the thermodynamic region (integer, or \texttt{None}) corresponding to the given temperature ($\degree$C) and pressure (Pa), as defined by the IFC-67 specification. The regions are:

\begin{enumerate}
  \item liquid water
  \item dry steam
  \item supercritical
  \item near-critical
\end{enumerate}

If the input temperature and/or pressure are outside the operating range of the IFC-67 formulation, the routine will return \texttt{None}.

\textbf{Parameters:}
\begin{itemize}
\item \textbf{t}: float\\
  Temperature ($\degree$C)
\item \textbf{Pressure}: float\\
  Pressure (Pa)
\end{itemize}

\section{Example}
\index{examples!thermodynamics}

The following script reads in a geometry file and writes an initial conditions file with approximate hydrostatic conditions corresponding to a specified vertical temperature gradient.  In this case, the model has a simple flat surface, so that each column has the same number of layers.  The \texttt{cowat} function is used to calculate the fluid density at each layer, and hence the approximate vertical pressure distribution.

\begin{lstlisting}
from mulgrids import *
from t2thermo import *

geo = mulgrid('gmodel.dat')

patm, tatm = 101.325e3, 15.0
ptblk = np.zeros((geo.num_blocks, 2))
ptblk[:,0] = patm; ptblk[:,1] = tatm

g = 9.8
p, t = patm, tatm
thick = 0.0
tgradient = 30 # deg C/km
for lay in geo.layerlist[1:]:
    d = cowat(t, p)[0]
    thisthick = lay.top - lay.bottom
    h = 0.5 * (thick + thisthick)
    p += d * g * h
    t += tgradient / 1.e3 * h
    thick = thisthick
    for col in geo.columnlist:
        blkname = geo.block_name(lay.name, col.name)
        iblk = geo.block_name_index[blkname]
        ptblk[iblk] = [p, t]
inc = dat.grid.incons(ptblk)
inc.write('model.incon')
\end{lstlisting}
